require(cjoint)
require(estimatr)

attr.names <- c("constitution", "attack", "pubwork", "constax", "corptax", 
                "surname", "samesex", "nuclear", "lockdown", "COVID")
level.names <- list(constitution = c("改正すべき", "改正すべきでない"), 
                    attack = c("検討すべき", "検討すべきでない"), 
                    pub.work = c("必要である", "必要でない"), 
                    cons.tax = c("引き下げるべき", "引き下げるべきでない"), 
                    corp.tax = c("引き上げるべき", "引き上げるべきでない"), 
                    surname = c("法律で認めるべき", "認めるべきでない"), 
                    samesex = c("法律で認めるべき", "認めるべきでない"), 
                    nuclear = c("廃止すべき", "維持すべき"), 
                    lockdown = c("法整備をすべき", "法整備をすべきでない"), 
                    COVID = c("治療体制の拡充", "治療薬やワクチンの開発・確保"))
attribute.labels <- c("憲法改正", "敵基地攻撃能力の保有", 
                      "公共事業による雇用確保", "消費税率", 
                      "企業が納める法人税率", "選択的夫婦別姓", 
                      "同性婚", "原子力発電", 
                      "強制的な外出自粛や休業措置", "コロナ対策の最優先事項")

#### データの読み込み・整形 ####
## 回答者レベルのデータ
respondent.data <- read.csv("main_data.csv")

## コンジョイントの候補者レベルのデータ
conjoint.data <- read.qualtrics("conjoint_data.csv", 
                                ranks = c("Q1_C_1", "Q1_C_2", "Q2_C_1", "Q2_C_2", 
                                          "Q3_C_1", "Q3_C_2", "Q4_C_1", "Q4_C_2", 
                                          "Q5_C_1", "Q5_C_2", "Q6_C_1", "Q6_C_2", 
                                          "Q7_C_1", "Q7_C_2", "Q8_C_1", "Q8_C_2", 
                                          "Q9_C_1", "Q9_C_2", "Q10_C_1", "Q10_C_2", 
                                          "Q11_C_1", "Q11_C_2", "Q12_C_1", "Q12_C_2", 
                                          "Q13_C_1", "Q13_C_2", "Q14_C_1", "Q14_C_2", 
                                          "Q15_C_1", "Q15_C_2"), 
                                covariates = "rid", new.format = TRUE)
# 変数名を英数字に
colnames(conjoint.data)[seq(6, 24, 2)] <- 
  c("COVID", "corptax", "lockdown", "constitution", "nuclear", 
    "pubwork", "constax", "surname", "samesex", "attack")
# コンジョイントの水準の変数を順序付き因子型に
for (i in 1:10) {
  conjoint.data[, attr.names[i]] <- factor(conjoint.data[, attr.names[i]], 
                                           level.names[[i]])
}

# 回答者レベルの比例投票先と長期的党派性の変数と結合
merged.data <- merge(conjoint.data, 
                     respondent.data[, c("PR", "PID", "rid")], 
                     by = "rid")
# 原稿提出後に再現コードを整理している時点で，この結合の処理は不要であることに
# 気付いたが，このプロセスを省いたコードを実行したところ，merge()を経る場合と
# 経ない場合とで，結論に影響するレベルではないものの，AMCEやIMCEの推定結果に
# ごくわずかな差が生じることが明らかになった。具体的に言えば，注23で報告した
# IMCEが有意になる回答者の割合の値が0.1ポイント異なるといったものである。
# このmerge()に関するRの挙動について，筆者には原因が思い当たらないが，掲載原稿
# に合わせてmerge()を使う過程を残しておくことにした。

#### AMCEの推定 ####
AMCE.est <- lm_robust(selected ~ constitution + attack + pubwork + 
                        constax + corptax + surname + samesex + 
                        nuclear + lockdown + COVID, 
                      data = merged.data, clusters = respondent)

# AMCEの推定結果（本文で言及）
round(summary(AMCE.est)$coefficients, 3)

# 「0.2ポイント程度」のAMCEの効果量（本文注21）
round(sd(merged.data$selected), 2)
round(-0.2 / sd(merged.data$selected), 2)
# なお，注21には本文で言及した「0.2ポイント程度」の効果量を記載したが，
# 実際の消費税率と同性婚のAMCEの効果量は次のとおりである。
round(summary(AMCE.est)$coefficients[5] / sd(merged.data$selected), 2)
round(summary(AMCE.est)$coefficients[8] / sd(merged.data$selected), 2)

# AMCEの推定結果の図（図A4）
cairo_pdf("Figure_A4.pdf", width = 4.5, height = 4, pointsize = 7, family = "Meiryo")
par(mar = c(3.5, 0, 0.5, 0.5), lwd = 0.5)
plot(NULL, NULL, type = "n", bty = "n", xlim = c(-0.7, 0.1), ylim = c(0, 30), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
abline(v = 0, col = "gray50")
abline(v = c(-0.3, -0.2, -0.1, 0.1), lty = 3, col = "gray50")
for (i in 1:10) {
  text(-0.7, 33 - 3 * i, attribute.labels[i], pos = 4)
  text(-0.7, 32 - 3 * i, 
       paste0("(比較対象: ", tolower(AMCE.est$xlevels[[i]][1]), ")"), pos = 4)
  segments(-0.31, 31 - 3 * i, 0.11, 31 - 3 * i, col = "gray80")
  text(-0.32, 31 - 3 * i, AMCE.est$xlevels[[i]][2], pos = 2)
  points(AMCE.est$coefficients[i + 1], 31 - 3 * i, pch = 19)
  segments(AMCE.est$conf.low[i + 1], 31 - 3 * i, 
           AMCE.est$conf.high[i + 1], 31 - 3 * i)
}
axis(1, at = seq(-0.3, 0.1, 0.1), labels = NA, lwd = 0.5)
mtext(c("-0.3", "-0.2", "-0.1", "0.0", "0.1"), 1, 1, 
      at = seq(-0.3, 0.1, 0.1))
mtext("AMCE", side = 1, at = -0.1, line = 2.5)
dev.off()

#### 属性間の交互作用 ####
# 交互作用項の係数の検定をする関数
interaction.t.test <- function(variable.1, variable.2, data) {
  lm.model <- lm_robust(as.formula(paste0("selected ~ constitution + attack + 
                                          pubwork + constax + corptax + surname + 
                                          samesex +nuclear + lockdown + COVID + ", 
                                          variable.1, " : ", variable.2)),
                        data = data)
  summary(lm.model)$coefficients[12, 4]
}

# i, j = 1, ..., 10, i > jとなる組み合わせについて検定
interaction.t.test.p.value <- matrix(NA, 10, 10)
for (i in 1:10) {
  for (j in 1:10) {
    if (i == j) break
    interaction.t.test.p.value[i, j] <-
      interaction.t.test(attr.names[i], attr.names[j],
                         data = merged.data)
  }
}

# 偽発見率補正後のp値（表A7）
round(matrix(p.adjust(interaction.t.test.p.value, method = "BH"), 10, 10), 3)
which(matrix(p.adjust(interaction.t.test.p.value, method = "BH"), 10, 10) < 0.05, 
      arr.ind = TRUE)

# 敵基地攻撃能力の保有と公共事業による雇用確保の交互作用の解釈
AMIE.est.1 <- lm_robust(selected ~ constitution + attack * pubwork + 
                          constax + corptax + surname + samesex + 
                          nuclear + lockdown + COVID, 
                        data = merged.data, clusters = respondent)
round(summary(AMIE.est.1)$coefficients, 3)
round(summary(AMIE.est.1)$coefficients[3] + 
        summary(AMIE.est.1)$coefficients[4] + 
        summary(AMIE.est.1)$coefficients[12], 3)

# 選択的夫婦別姓と同性婚の交互作用の解釈
AMIE.est.2 <- lm_robust(selected ~ constitution + attack + pubwork + 
                          constax + corptax + surname * samesex + 
                          nuclear + lockdown + COVID, 
                        data = merged.data, clusters = respondent)
round(summary(AMIE.est.2)$coefficients, 3)
round(summary(AMIE.est.2)$coefficients[7] + 
        summary(AMIE.est.2)$coefficients[8] + 
        summary(AMIE.est.2)$coefficients[12], 3)

#### IMCEの推定 ####
## IMCEの推定
N <- max(merged.data$respondent)
IMCE.result <- data.frame(matrix(NA, N, 66))
colnames(IMCE.result) <- c("rid", 
                           paste0(c("intercept", attr.names), ".est"), 
                           paste0(c("intercept", attr.names), ".se"), 
                           paste0(c("intercept", attr.names), ".lower"), 
                           paste0(c("intercept", attr.names), ".upper"), 
                           paste0(c("intercept", attr.names), ".p"), 
                           paste0(attr.names, ".padj"))
for (i in 1:N) {
  # 回答者個人のデータを取り出す
  individual.data <- subset(merged.data, respondent == i)
  # IMCEを推定
  IMCE.est <- lm_robust(selected ~ constitution + attack + pubwork + 
                          constax + corptax + surname + samesex + 
                          nuclear + lockdown + COVID, 
                        data = individual.data)
  # 回答者ID
  IMCE.result$rid[i] <- individual.data$rid[1]
  # 点推定値
  IMCE.result[i, 2:12] <- IMCE.est$coefficients
  # 標準誤差
  IMCE.result[i, 13:23] <- IMCE.est$std.error
  # 95%信頼区間
  IMCE.result[i, 24:34] <- IMCE.est$conf.low
  IMCE.result[i, 35:45] <- IMCE.est$conf.high
  # p値
  IMCE.result[i, 46:56] <- IMCE.est$p.value
}

# 偽発見率によるp値の補正
for (i in 1:10) {
  IMCE.result[, 56 + i] <- p.adjust(IMCE.result[, 46 + i], "BH")
}

## IMCEの推定値の分布（図2）
cairo_pdf("Figure_2.pdf", width = 6, height = 5, pointsize = 7, family = "Meiryo")
layout(matrix(1:12, 3, 4, byrow = TRUE))
par(mar = c(0.5, 0, 3, 0), lwd = 0.5)
for (i in 1:10) {
  plot(NULL, NULL, type = "n", bty = "n", xlim = c(-5.7, 5.7), ylim = c(-3, 1.2), 
       xlab = "", ylab = "", xaxt = "n", yaxt = "n")
  segments(-5:5, rep(0, 11), -5:5, rep(-3, 11), lty = 3, col = "gray")
  segments(-5, 0, 5, 0)
  points(-5:5, rep(0, 11), pch = "|", cex = 0.7)
  text(c(-5, 0, 5), rep(0, 3), c("-5", "0", "5"), pos = 3)
  density.est <- density(IMCE.result[, i + 2])
  lines(density.est$x, density.est$y)
  temp.data <- IMCE.result[order(IMCE.result[, i + 2]), 
                           c(i + 2, i + 24, i + 35, i + 46)]
  y <- seq(-3, -0.1, length.out = nrow(temp.data))
  for (j in 1:nrow(temp.data)) {
    segments(temp.data[j, 2], y[j], temp.data[j, 3], y[j], lwd = 0.01, 
             col = ifelse(temp.data[j, 4] < 0.05, "gray30", "gray80"))
  }
  lines(temp.data[, 1], y)
  title(paste0(attribute.labels[i], "\n(", level.names[[i]][2], ")"), 
        line = 0.5, cex.main = 1.2)
}
dev.off()

## p値が0.05を下回った人の割合（本文で言及）
for (i in 1:10) {
  cat(attr.names[i], ":\n")
  print(round(mean(IMCE.result[, paste0(attr.names[i], ".p")] < 0.05), 3))
}

## 補正後のp値が0.05を下回った人の割合（注23，24）
for (i in 1:10) {
  cat(attr.names[i], ":\n")
  print(round(mean(IMCE.result[, paste0(attr.names[i], ".padj")] < 0.05), 3))
}

#### IMCEと投票政党・党派性の関係 ####
## IMCEの推定結果と回答者レベルのデータを結合
IMCE.data <- merge(respondent.data, IMCE.result, by = "rid")

# 小政党をその他にまとめる
IMCE.data$PR.2 <- ifelse(IMCE.data$PR < 8, IMCE.data$PR, 
                         ifelse(IMCE.data$PR == 11, 8, 9))
IMCE.data$PID.2 <- ifelse(IMCE.data$PID < 8, IMCE.data$PID, 
                          ifelse(IMCE.data$PID == 11, 8, 9))

## 比例投票先別のIMCEの推定値の分布（図A5）
PR.labels <- c("自民", "立憲", "公明", "共産", 
               "維新", "国民", "れいわ", "棄権")
for (i in 1:10) {
  cairo_pdf(paste0("Figure_A5_", i, ".pdf"),
            width = 8, height = 0.9, pointsize = 6, family = "Meiryo")
  layout(matrix(1:8, 1, 8))
  par(mar = c(2, 0.2, 0, 0.2), oma = c(0, 0, 3, 0), lwd = 0.5)
  for (j in 1:8) {
    plot(NULL, NULL, type = "n", bty = "n", xlim = c(-2, 2), ylim = c(0, 1.5), 
         xlab = "", ylab = "", xaxt = "n", yaxt = "n")
    density.est.1 <- density(IMCE.data[IMCE.data$PR.2 == j, 
                                       paste0(attr.names[i], ".est")])
    density.est.2 <- density(IMCE.data[IMCE.data$PR.2 != j, 
                                       paste0(attr.names[i], ".est")])
    polygon(density.est.2 $x, density.est.2 $y, border = NA, col = "gray80")
    lines(density.est.1$x, density.est.1$y)
    if (i %% 5 == 1) {
      text(-2, 1.4, PR.labels[j], pos = 4, cex = 1.5)
    }
    axis(1, lwd = 0.5)
  }
  title(paste0(attribute.labels[i], " (", level.names[[i]][2], ")"), 
        cex.main = 1.5, outer = TRUE, line = 1)
  dev.off()
}

## 長期的党派性別のIMCEの推定値の分布（図A6）
PID.labels <- c("自民", "立憲", "公明", "共産", 
                "維新", "国民", "れいわ", "無党派")
for (i in 1:10) {
  cairo_pdf(paste0("Figure_A6_", i, ".pdf"),
            width = 8, height = 0.9, pointsize = 6, family = "Meiryo")
  layout(matrix(1:8, 1, 8))
  par(mar = c(2, 0.2, 0, 0.2), oma = c(0, 0, 3, 0), lwd = 0.5)
  for (j in 1:8) {
    plot(NULL, NULL, type = "n", bty = "n", xlim = c(-2, 2), ylim = c(0, 1.5), 
         xlab = "", ylab = "", xaxt = "n", yaxt = "n")
    density.est.1 <- density(IMCE.data[IMCE.data$PID.2 == j, 
                                       paste0(attr.names[i], ".est")])
    density.est.2 <- density(IMCE.data[IMCE.data$PID.2 != j, 
                                       paste0(attr.names[i], ".est")])
    polygon(density.est.2 $x, density.est.2 $y, border = NA, col = "gray80")
    lines(density.est.1$x, density.est.1$y)
    if (i %% 5 == 1) {
      text(-2, 1.4, PID.labels[j], pos = 4, cex = 1.5)
    }
    axis(1, lwd = 0.5)
  }
  title(paste0(attribute.labels[i], " (", level.names[[i]][2], ")"), 
        cex.main = 1.5, outer = TRUE, line = 1)
  dev.off()
}